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Interpolation is often used to improve the 
accuracy of integrals over spectral data 
convolved with various response functions 
or power distributions. Formulae are 
developed for propagation of uncertainties 
through the interpolation process, specifi- 
cally for Lagrangian interpolation increas- 
ing a regular data set by factors of 5 and 2, 
and for cubic-spline interpolation. The 
interpolated data are correlated; these cor- 
relations must be considered when com- 
bining the interpolated values, as in inte- 
gration. Examples are given using a com- 
mon spectral integral in photometry. 



Correlation coefficients are developed for 
Lagrangian interpolation where the input 
data are uncorrelated. It is demonstrated 
that in practical cases, uncertainties for the 
integral formed using interpolated data can 
be reliably estimated using the original 
data. 
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1. Introduction 

Measurements of spectral irradiance, spectral 
responsivity or spectral reflectance are often made for a 
limited set of wavelengths and then used to calculate 
weighted spectral sums for photometry, colorimetry or 
filter radiometry. It is often necessary to interpolate the 
spectral data to a finer grid to avoid errors arising from 
the discrete approximation used to estimate the integral 
where the weighting function varies strongly between 
the wavelengths at which the measurements were made 
[1]. Interpolated values are then correlated to the origi- 
nal data and to nearby interpolated points; unless these 
correlations are taken into account, uncertainty calcula- 
tions will give misleading results, generally under- 
estimating errors in the spectral sums by significant 
amounts. 

Interpolation may also be required when reference 
standards are provided for a limited set of wavelengths. 
Many of the primary reference standards in the nation- 
al metrology institutes are derived in a functional form, 
and can be calculated on as fine a wavelength grid as 
required, along with all the necessary correlations [2]. 
However, calibrations of client lamps or detectors 
involve a transfer, by comparison, to artifacts that do 



not necessarily show a spectral variation that is easily 
modeled, and in the interests of reducing costs may be 
provided on a limited number of wavelengths. The 
reference data at different wavelengths may also be 
correlated. While the primary reference standards are 
often strongly correlated between wavelengths, the 
transfer process itself adds uncertainty that is generally 
random and often reduces the correlations to negligible 
levels [2] 

Where data are available at sufficient wavelengths to 
avoid errors due to the sum approximation to the inte- 
gral, it is preferable to interpolate reference tables, such 
as the photometric response function V x and the colori- 
metric tristumulus response functions [3] to the wave- 
lengths at which measurements are taken, because 
those tables contain no uncertainty However, inter- 
polation of measurement data is often applied. One 
reason is that software for a calculation may require 
data at a particular interval, another is that instrument 
measurement programs may provide limited sets of 
data. 

Following a brief description of uncertainty propaga- 
tion, this paper is divided into two sections. The first 
covers Lagrangian interpolation, the second cubic 
spline interpolation; these are the two most commonly 
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used interpolation methods. In each section, interpola- 
tion from a data set which is itself correlated is consid- 
ered. Various simplifications for practical applications 
are made, and examples are presented. A conclusion is 
that in practical terms, uncertainty can be accurately 
derived from the original data set without a complex 
calculation of correlations. 

2. Propagation of Uncertainty 

Uncertainty propagation is described in detail in the 
ISO Guide to the Expression of Uncertainty in 
Measurement [4]. The uncertainty in a quantity y 
formed by combining measured quantities x { through 
the relationship y =f(x h x 2 , ..%) is given by 



u\y): 



if— 

1^ dx, 



» 2 w+££^^«(**,)> a) 



=i jj4=i dx t uXj 



where wfo) is the uncertainty in jc,- and u(x h xj) is the 
covariance between x t and x r For uncorrelated input 
quantities, the covariance between pairs of variables is 
zero and Eq, (1) reduces to the "sum of squares" com- 
monly applied. The derivatives 3/7 3x, are sensitivity 
coefficients for the dependence of y on the various 
measured quantities. Given that u 1 (x,) = u(x h x t ), 
Eq. (1) can be expressed as 



u\y) = f y U x f i „ 



(2) 



where 



/= 



3/3/ 



Y 



(3) 



3x 3x 9 dx, 

y 1 2 ; 

is a column vector of sensitivity coefficients ( T indi- 
cates the transpose) and 

U x = [u( Xi , Xj )] (4) 

is the Nx TV uncertainty matrix. 

Interpolation of spectral data is generally performed 
to produce a set of values (p k ,x k ) from set (y t , x t ) 
where x t is the independent variable (generally wave- 
length in radiometry). Quantities in the set (p k ,x k ) that 
depend on the same (y t , x t ) are correlated through this 
common dependence. The covariance between two 
values p k and p m (and, when k = m, the square of the 
uncertainty) is given by 



u{p k ,P m )=tt^ d f-u{ yi , yj ). 
Mjioy, tyj 



(5) 



In matrix form, this is simply expressed as 

u{P k ,P m ) = f Pk T U y f„„- (6) 

It is sometimes convenient to use correlation coeffi- 
cients rather than covariances, defined as 

<Pk,P m ) 



r{Pk,P m ) = 



J»(Pt,PtMPm>Pj 



(7) 



A matrix of correlation coefficients is square, symmet- 
ric about the diagonal and has value 1 in the diagonal 
elements. 

3, Lagrangian Interpolation 

We have a tabulated function y i at values x;. 

yi =y{ Xi ),i = \....N. (8) 

For x in the range x m to x m+n , the formula for 
Lagrangian interpolation is [5] 

n 

. n (*-*«.*- 1) 



1 1 \ X m+j-l X m+i-\ ) 



(9) 



: £^j m+ ^ 



7=1 



Equation (9) represents an (w-l)th order polynomial fit- 
ted through the n original values. Interpolated data are 
formed as a linear combination of nearby existing data. 
Sensitivity coefficients for the dependence of the inter- 
polated data on the input data are simply the weights Wj 
in Eq. (9). Covariances between the output values, that 
is, the original input values plus the interpolated values, 
take several forms. Correlations present in the original 
values remain. The values newly formed by interpola- 
tion are correlated to both the input values forming 
them (and through them to the remaining input values 
if correlations are present in the original set) and to any 
new values formed from common input values. All of 
these covariances, including the uncertainty of the 
interpolated values, can be calculated through Eq. (6). 
In many instances, the input data y- t are present at 
regular values of x t and further, the output data are also 
required at regular intervals. Interpolation is usually 
then performed by running Eq. (9) through the data and 
forming a new value or values in the center of the 
range. The multipliers convolving the data are deter- 
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mined by the interval spacings only and can be calcu- 
lated prior to the interpolation. This was clearly demon- 
strated by Savitsky and Golay [6] in the development of 
their algorithms for smoothing and differentiation of 
spectral data where polynomial expressions of various 
order are fitted to regularly- spaced data; in those cases, 
the coefficients are determined as fixed linear combina- 
tions of the input dependent data. The arguments pre- 
sented here for Lagrangian interpolation can easily be 
extended to cover smoothing of data using the 
Savitsky-Golay routines. Two-point Lagrangian inter- 
polation forming a new value in the center of the exist- 
ing values has weights (w h w 2 )=( Vi, V2), equivalent to a 
linear interpolation. 

Propagation of uncertainty through two common 
examples of Lagrangian interpolation used in photo- 
metry are now discussed. In both of these we consider 
the calculation of illuminance response to CIE illumi- 
nant D65, a tabulated distribution carrying no uncer- 
tainty, for a photometer whose spectral response is a 
close approximation to V A , measured at different 
spectral intervals and where the measured response 
values are uncorrected. These distributions are shown 
in Fig, 1 for a wavelength interval of 5 nm. As the 
response function tapers to zero at each end, the 
luminance response is given as 



R v -AA2^£ D65f ., 



(10) 



where R ( is the photometer response and E D65J is the 
illuminant value at the zth wavelength, respectively, and 
AX is the 5 nm wavelength separation between the 
values. 



For uncorrelated spectral response values, the uncer- 
tainty in R v is given by 

u\R r ) = (AX) 2 f j E 2 D65 .u\B i ). (11) 



For values as tabulated by CIE [3], the value of R v is 
10567.41, with an uncertainty 19.60 (relative uncer- 
tainty 0.1855 %) if the responsivity values have a rela- 
tive uncertainty of 1 % and are uncorrelated. (Note that 
extra significant figures for the uncertainty are present- 
ed above those of normal practice for the purpose of 
comparison.) 

3.1 Photometer Measured at 10 nm 
Intervals, Interpolated to 5 nm 

The spectral integral Eq. (10) for input values on a 
10 nm grid evaluates to 10568.18, a small change com- 
pared to the 5 nm data due to the discrete approxima- 
tion to the integral; the relative uncertainty for uncorre- 
lated spectral response values with a relative uncertain- 
ty of 1 % becomes 0.2623 %, or the expected V2 
increase compared to the response measured at 5 nm 
intervals. We wish to interpolate the photometer 
response with a four point Lagrangian function to data 
on a 5 nm interval. The weights for Eq. (9) are then 



W l9 W 2 ,Ws,W 4 



-1 



-1 



16 16 16 16 



(12) 



and, as the input data are uncorrelated, the uncertainty 
for an interpolated value is given by 
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Fig. 1. Distributions CIE V^ and D65 used for interpolation examples. 



In any given interval spanning four input values, the 
uncertainty value u (y 7 ) is approximately constant, and 
the interpolated values have an uncertainty approx. 80 
% of the original values in that range. If we ignore the 
correlations that have been introduced, the relative 
uncertainty of the integral evaluates as 0.17%, low 
compared to the original value and clearly incorrect as 
it would imply that it could be reduced to zero by 
repeated interpolation. 

A full correlation matrix for the output values is 
formed as follows. First the interpolation is performed 



71 



Volume 108, Number 1, January-February 2003 

Journal of Research of the National Institute of Standards and Technology 



(14) 



(15) 



(using linear interpolation in the first and last intervals 
of the input values). From TV original values, we now 
have (27V- 1) values 
3'/=3'oriBhud,^'odd J 7=(i+l)/2 

yi = ^1^-3+^2^-1 +w 3 # +1 w 4 # +3 ; 

4<i<2N-4, /even. 
Uncertainties for these values are known for the origi- 
nal data and for the interpolated data from Eq. (13) and 
these are used to populate the diagonal elements of the 
uncertainty matrix U y We then have to populate the 
elements to the right of the diagonal only before filling 
to the left of the diagonal by symmetry. As the original 
input values (now at / odd) are uncorrected, we have 
for all /, from Eq. (5), 

«0 ; /»3 ; / + i) = W3tt 2 U +1 ) 

m (^i^i +3 ) = w 4 m2 (^ + 3) 
For / even (interpolated values), 

y i+ i = w i^_i + w 2 yi + i + w 3^ +3 + ^y i+5 

y i+ 4 = ^+1+^+3 +^^5+^4^7 > 

y i+ 6 = w iyt + z +w 2yM + ^y i+7 +^ 4 ^ 

and from Eq, (5) we have 

u (yi>y i+ 2) = "i^CX-i) + ™> w*u 2 (y i+1 ) + w 3 w 4 u(y i+3 ) 

u(y i9 y i+A ) = ^i^u{y i+l ) + w 2 w 4 u 2 (y l+3 ) 

u(y i9 y i+6 ) = w 1 w 4 u 2 (y i+3 ) 

(17) 

For completeness, similar expressions were applied for 
the values formed by linear interpolation in the first and 
last intervals. Uncertainty calculated for the integral 
including all these correlations was then exactly that 
calculated with the original data for the 10 nm grid. 

One further simplification can be made for uncorre- 
cted input values. In practical terms, nearby values 
have the same uncertainty. Hence the sets of Eqs. (14) 
to (17) can be reduced to a matrix of correlation coeffi- 
cients determined purely from the weights, 
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where the first row corresponds to an interpolated 
value. Uncertainty for the integral using the interpolat- 
ed data set is then found by modifying the sensitivity 
column vector to include the uncertainty at each value, 



/= 






«Oi) 



(19) 



and then performing the matrix multiplication Eq. (2). 
A negligible change relative to the true values is due to 
averaging through regions where the response is chang- 
ing rapidly; while the relative uncertainty in these 
regions is constant, the absolute value is not. 

Table 1 shows uncertainties calculated for all the 
options discussed in this section. It can be seen that 
proper accounting for the correlations introduced by the 
interpolation reproduces the uncertainty calculated for 
the input values alone, and that using the correlation 
coefficients provides a practical calculation of the 
uncertainty matrix. 

Table 1. Integral value and relative uncertainty for various calcula- 
tion options; four point Lagrange interpolation adding one value in 
each range. Input values uncorrelated, 1 % relative uncertainty. 



Method 


Integral 


Relative uncertainty 


Original data on 5 nm grid 


10567.4 


0.185 


Original data on 10 nm grid 


10568.2 


0.262 


Interpolated data; ignore 


10567.5 


0.168 


correlations 






Interpolated data; with 


10567.5 


0.262 


correlations 






Interpolated data; use 


10567.5 


0.266 


correlation coefficients 







3.2 Photometer Measured at 5 nm 
Intervals, Interpolated to 1 nm 

After re-arranging the N input values to their new 
positions, we add four values between input values y, 
and y l+5 that can be represented as 

r -6 108 27 -4 v 
1 -8 84 56 -7 
125 -7 56 84 -8 
-4 27 108 -6 



'y M ^ 



y i+2 
y M 
y + 4 



'y* " 



y 

y + s 

y + m 



■ (20) 



Fig. 2 shows the relative uncertainty of the central data 
set where a photometer with a V x response is interpolat- 
ed from 5 nm to 1 nm and the input values are assumed 
uncorrelated with a relative uncertainty of 1 %. 
Uncertainties of the interpolated values are lower than 
those of the input values. If we ignore the correlations 
introduced by the interpolation, the relative uncertainty 
of the integral with the D65 illuminant (interpolated to 
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a 1 nm grid with a cubic-spine routine), reduces to 
0.074 %, a value too low by the order of Vs . The 
integral itself has the same value as that shown in 
Table 1 for the 5 nm grid. 

For a four point Lagrange interpolation adding four 
values between each of the input values, correlations in 
the output set extend over the 19 values following an 
input value. The matrix of correlation coefficients is 
shown as the transpose relative to Eq. (18) in the inter- 
ests of printing; that is, it is equivalent to filling to the 
lower left of the diagonal of the correlation matrix prior 
to filling the upper right by symmetry. Beginning at a 
column corresponding to an input value, the matrix of 
correlation coefficients is 
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Fig. 2. Relative uncertainties for the interpolated V^ data set. 
Original values assumed uncorrected, with a relative uncertainty of 
1%. 
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A number of these, for the values furthest from the 
input set, are negligible. The relative uncertainty of the 
integral for the interpolated data set calculated with 
these correlation coefficients is 0.197 %, equivalent in 
practical terms to that calculated using the original 
5 nm data set and shown in Table 1 . 

4. Cubic-Spline Uncertainty Propagation 

For our set of data Eq. (8), cubic-spline interpolation 
[7] calculates a value y at x in the interval x t to x m as 



where 



y = Ay i +By i+1 + Cy i +Dy l+1 , 



A=^±^, 



B- 



(22) 

(23) 
(24) 



Then even for uncorrected input data, the output data 
are correlated over the whole set of interpolated values. 
We wish to calculate the co variance u (y n , y m ) between 
two interpolated values 



y m = Ay* + B m y i+ i + C J\ + D m i + i 
y n = Ayj + B n y J+1 +c n y] + Dj j+l 



> (29) 



where the y n and y m values may be in the same or differ- 
ent intervals denoted by ij. The uncertainty in y n is 
given by 



u 2 (y n ) = u(y H ,y H ). 



(30) 



The covariance between (and uncertainty of) the input 
values is known, carried in the matrix 



U,=[u(y t , yj j\. 



(31) 



C = -(A 3 -A)(x M - Xl ) 2 , 
6 

D=±(B 3 -B)(x i+1 - Xi ) 2 . 
6 



(25) 
(26) 



The first two terms of Eq. (22) represent simple linear 
interpolation. Including the second derivatives y" 
yields a function that has first and second derivatives 
continuous at the boundaries between intervals. 

The second-derivatives are unknown. The relation 
between them is given by 



X- X- i « X-,-\ X- i " X-.-i X- « 

6 3 6 



yj+i yt yi y^- 



(27) 



X i+l X i X i X i-\ 



which is a system ofN-2 equations in the A^ unknowns 
y. . The natural cubic-spline, which is commonly used, 
sets 



yt=y m N=Q 



(28) 



and solves for the remaining terms [7]. We are interest- 
ed in using the cubic-spline interpolation on spectral 
data of known uncertainties, including the possibility of 
correlations, where the interpolated data may then be 
combined in various ways, so that not only the uncer- 
tainties in the interpolated data but also the correlations 
present are important in propagating uncertainties in 
the combinations. 

The (iV-2) values of y t depend on each of the 
input values, i.e., are correlated to each input value. 



To propagate uncertainties through Eq. (29) we need 
the covariance between the second-derivative values 
and the input values, and between the second-deriva- 
tives themselves. These in turn require the sensitivity 
coefficients for the dependence of the second-deriva- 
tive on the input values, dy k ldy t . 

The set ofN-2 equations Eq. (27) can be written in 
matrix form 



J3-J2 y 2 -yi 



y i+ 2-y i+ i y M -yt 



y N -y N -i y N -i-y N -2 



= (Kj) 



y 2 



y N - 



(32) 



where 



*JV-1 **JV-2 



^ =iZ f ii 



h, M = 



6 

X i+2 ~ X i 

3 

X; , ~, X;, 



(33) 



_ -*i+2 ^/+1 



is a tri-diagonal matrix with remaining terms zero. 
Multiplying both sides of Eq. (32) on the left by 

the inverse matrix (\ 7 ) yields the required second- 
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derivatives. By selecting the Ufa row of this multiplica- 
tion we have the relationship 



y* 



■IK 



f y J+ 2-y J+ x y j+ i-yj ^ 



7=1 



(34) 



which when rearranged to 



y i+ 



7=3 

N-l ( 

7=2 

N-2 



y j+ i 



+ 



y } 



(35 



j=i 



yields 



Ma 



KU 



-h. 



ij-l 



K, 



~j+i 



These are augmented by 

dy k /dy i =0,k = l,N. 



(36 



(37) 



If the input y l values are uncorrected, the required 
covariances are given by 






(38) 



and 



where 



u(y;, yj ) = ±f^u\y t ) 

"mm 






u 2 (y k ) 



(39) 



Note that the NxN matrix represented by Eq. (39) is 
not symmetric. For correlated inputs, the sums required 
are the matrix products 



»(y i ,y J ) = f?u y f J 
»(y~ t ,yj) = f?u ygj 



(40) 



where 



fr 



3y* 



,k = \..N 



(41) 



is a column-vector of sensitivity coefficients of the 
second-derivatives vs the input values and gj is a 
column vector of length N with 1 in the y'th row, 
elsewhere. 

Covariances and uncertainties for the interpolated 
values are then found by recognizing that the values of 
y m and>> tt in Eq. (29) are represented by 



" 



f m {yi y M y] y M yj y J+ i y) y J+ if ( 42 ) 



where 



Jn j n yyi -^I'+i y% yi+i yj yj+i yj yj+i* 



f m =(A i B i C i D i 0) 
/.=(0 OAjBjCjDj) 



(43) 



are sensitivity vectors for y m and y„ against the eight 
variables in the vectors shown in Eq. (42). Covariances 
between the interpolated values are then given by 



<y m ,y n ) = f m Uf n9 



(44) 



U = 



r u{y t ) u(y i+l ,y t ) u(y,y) u(y M9 y) u{y p y) i(y j+19 y) u(y j9 y) t(y J+1 , y) 
u{y i+l ) u(y n y M ) u{y i+l9 y i+l ) u{y p y i+l ) *0, +1 , y i+1 ) u(y p y +1) t^y J+1 , y l+1 ) 

u 2 (y°) u{y i+l ,y t ) u{y p y) u(y j+v y) u(y p y) u{y j+l ,y t ) 

u 2 {y i+l ) u{y p y i+l ) u(y j+l , y i+1 ) u{y p y i+1 ) u(y J+1 , y i+1 ) 

uiyj) u(y j+19 yj) u{y p yj) u(y j+19 yj) 

u 2 {y j+l ) <y p y j+l ) u{y j+l9 y J+1 ) 

u 2 {y)) u(y] + i>y]) 
u 2 {y j+x ) 



(45) 



75 



Volume 1 08, Number 1 , January-February 2003 

Journal of Research of the National Institute of Standards and Technology 



with the lower half symmetric about the diagonal. Only 
half of this matrix is required. The two right quadrants 
can be separately multiplied on the right by the column 
vector (Aj Bj C } D y ) T and the two results combined 
into a single eight-element column vector for the final 
multiplication. 

4.1 Cubic-Spline Interpolation Examples 

Consider again the photometer response curve of 
Fig. 1 , to be integrated over the wavelength range from 
360 nm to 830 nm (effectively the photometric 
response to an equal-energy source). The function was 
interpolated over the same input range (one less value) 
but shifted by 2.5 nm and the integral recalculated. 
Similarly the function was interpolated to 1 nm inter- 
vals and the integral recalculated. Table 2 shows the 
results for these interpolations, where the uncertainty in 
the integral was calculated based on 1 % uncertainty in 
the input values, uncorrected between values. The con- 
sequence of ignoring the correlation between the inter- 
polated values is also shown. Correlations between dis- 
tant points, introduced through the dependence of the 
second-derivatives, were negligible (largely because 
the response curve is relatively smooth), but strong cor- 
relations were found between near-neighbours. 

Table 2. Integral of V x and its uncertainty for input values on a 5 nm 
grid and 1 % relative uncertainty, uncorrected 



Figure 3 shows propagated uncertainties for the inter- 
polation shifting the input by 2.5 nm; for an interpola- 
tion to the mid-point, we would expect the 
interpolated value for a smooth function to be near the 
mean of the two interval boundaries with a propagated 
uncertainty 1/72=0.71 of the input (but of course 
correlated to adjacent values). 

Figure 4 shows the variation in uncertainty for values 
interpolated at different positions within the 5 nm inter- 
val. This is a practical concern where a wavelength off- 
set may be present in measurements, although in gener- 
al it is a better practice to retain the wavelength values 
of the measured points and interpolate weighting func- 
tions such as the illuminant or the colorimetric response 



functions as these carry no uncertainty. At the input 
points, uncertainties equal that of the input; at the mid- 
points they fall tol/V2 of the input points. Again this 
is expected as the input data are smooth, and cubic- 
spline results are not much different from linear inter- 
polation. The cubic-spline reproduces the input ordi- 
nate values for abscissa values equal to an input value; 
for these points, we expect the propagated uncertainty 
to be that of the input and correlations between similar 
such values to match that of the input matrix U y . These 
conditions can be used to test the coding. 

Figure 5 shows V k interpolated from a 20 nm grid to 
2 nm. Where the input curve is changing rapidly rela- 
tive to the magnitude of the data, linear interpolation 
would be discontinuous and for these regions, the rela- 
tive uncertainty of the interpolated values rises above 
that of the 1 % assumed for the input values. This is 
shown more strongly in Fig. 6 for input data spaced at 
40 nm, interpolated to 5 nm. Here the interpolation 
does not provide a good representation of Fx(as shown 
in Fig. 7); where the interpolation is poor, the uncer- 
tainties for the interpolated values rise above those of 
the input. 
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5, Conclusion 



Interpolation of spectral data is a common occur- 
rence in radiometric and photometric measurements. 
Those data are often then combined in forming integral 
values such as photometric or colorimetric responses or 
filter radiometer responses. Interpolation is particularly 
important when a relatively smooth curve available 
only on a wide spectral spacing may need to be con- 
volved with a more-rapidly changing curve and then 
integrated. Uncertainties in the interpolated values will 
generally be smaller than those of the input data, 
although this is not always true in the case of cubic- 



spline interpolation. Uncertainties in combinations 
calculated using the interpolated data set then must 
include correlations introduced by the interpolation . 
Ignoring the correlation will lead to significant under- 
estimation of uncertainties. The calculations in this 
paper for both Lagrange interpolation and for cubic- 
spline interpolation show that the uncertainty can be 
reliably estimated, in practical terms, by propagating 
the uncertainty through the combination using only the 
original set of data. 
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